# Librarieslibrary(tidyverse)library(gmRi)library(rnaturalearth)library(scales)library(sf)library(gt)library(patchwork)library(lmerTest)library(emmeans)library(merTools)library(tidyquant)library(ggeffects)library(performance)library(gtsummary)library(gt)library(sizeSpectra)library(ggdist)library(pander)library(ggh4x)library(nlme)library(ggimage)conflicted::conflict_prefer("select", "dplyr")conflicted::conflict_prefer("filter", "dplyr")conflicted::conflicts_prefer(lmerTest::lmer)# Processing functionssource(here::here("Code/R/Processing_Functions.R"))# ggplot themetheme_set(theme_gmri(axis.line.y =element_line(),axis.ticks.y =element_line(), rect =element_rect(fill ="white", color =NA),panel.grid.major.y =element_blank(),strip.text.y =element_text(angle =0),axis.text.x =element_text(size =8),axis.text.y =element_text(size =8),legend.position ="bottom"))# labels for factor levelsarea_levels <-c("GoM", "GB", "SNE", "MAB")area_levels_long <-c("Gulf of Maine", "Georges Bank", "Southern New England", "Mid-Atlantic Bight")area_levels_all <-c("Northeast Shelf", area_levels)area_levels_long_all <-c("Northeast Shelf", area_levels_long)# table to join for swapping shorthand for long-hand namesarea_df <-data.frame(area =c("Scotian Shelf", "Gulf of Maine", "Georges Bank", "Southern New England", "Mid-Atlantic Bight", "Northeast Shelf"),survey_area =c("SS", "GoM", "GB", "SNE", "MAB", "Northeast Shelf"),area_titles =c("Scotian Shelf", "Gulf of Maine", "Georges Bank", "Southern New England", "Mid-Atlantic Bight", "Northeast Shelf"))# Degree symboldeg_c <-"\u00b0C"# Function to get the Predictions# Flag effect fits that are non-significant ###get_preds_and_trendsignif <-function(mod_x){ modx_preds <-as.data.frame(# Model predictionsggpredict( mod_x, terms =~ est_year * area * season) ) %>%mutate(area =factor(group, levels = area_levels_long_all),season =factor(facet, levels =c("Spring", "Fall")))# Just survey area and year modx_emtrend <-emtrends(object = mod_x,specs =~ area*season,var ="est_year") %>%as_tibble() %>%mutate(zero =0,non_zero =if_else(between(zero, lower.CL, upper.CL), F, T))# Preds with signif modx_preds %>%left_join(select(modx_emtrend, area, season, non_zero))}
Storycrafting - Northeast Shelf Spectra
I think we need to step back a bit and find the main story again–which is about how size structure of the NES fish community has changed over time, whether that varies by region, and what covariates are associated with those changes.
The Aug draft you developed really focused on comparing regression models to consider drivers based on the fish community size spectrum as viewed from the perspective of the Wigley species set vs. the full species set.
However, I didn’t see figures of those spectra in that draft (maybe there are somewhere else and I’ve lost track). The recent slide deck really focuses on seasonal differences in the size spectra and influence of temperature. - KMills
Everything in this markdown uses the size spectra exponent values where I shifted the wmin parameter to 16g. The wmin values are fixed at 16g across all groups.
Size spectra are steepening in 3 seasonal communities. Gulf of Maine (Spring), Georges Bank (Spring), and Southern New England (Spring).
The size distribution in the Mid-Atlantic Bight (Spring) is becoming more shallow.
Code
# Model linear trends in timespectra_trend_mod <-gls( b ~scale(est_year) * area * season,correlation =corAR1(form =~ est_year | area/season),data = wigley_bmspectra_model_df )#check_model(spectra_trend_mod)# plot(check_collinearity(spectra_trend_mod))# Get the predictions and flag whether they are significantbmspectra_trend_predictions <-get_preds_and_trendsignif(spectra_trend_mod) %>%mutate(model_id ="bodymass_spectra-Wigley Subset") %>%separate(model_id, into =c("metric", "community"), sep ="-") %>%mutate(metric =fct_rev(metric),area =factor(area, levels = area_levels_long_all))# Make the plotbmspectra_trend_p <-ggplot() +geom_ribbon(data =filter(bmspectra_trend_predictions, non_zero ==TRUE),aes(x, ymin = conf.low, ymax = conf.high, fill = season),linewidth =1, alpha =0.35) +geom_point(data = wigley_bmspectra_model_df ,aes(est_year, b, color = season),size =0.8, alpha =0.5) +geom_ma(data = wigley_bmspectra_model_df ,aes(est_year, b, color = season, linetype ="5-Year Moving Average"), n =5, ma_fun = SMA, key_glyph ="timeseries" ) +geom_line(data =filter(bmspectra_trend_predictions, non_zero ==TRUE),aes(x, predicted, color = season, linetype ="Significant Trend"),linewidth =1, alpha =1) +scale_fill_gmri() +scale_color_gmri() +facet_grid(area ~ ., scales ="free") +labs(x ="Year",y ="Size Distribution Exponent",linetype ="Trend",color ="Season",fill ="Season")# Show plotbmspectra_trend_p
Median body weight is declining in 3 seasonal communities, Gulf of Maine (Spring and Fall) and Georges Bank (Spring).
Median weight in the Mid-Atlantic Bight has been increasing in the Spring.
Code
# Load the median weight/length datawigley_size_df <-read_csv(here::here("Data/processed/wigley_species_size_summary.csv"))# shelf-wide summarywigley_sizes_shelf <-read_csv(here::here("Data/processed/shelfwide_wigley_species_size_summary.csv"))# Join themsize_results <-bind_rows(wigley_size_df, wigley_sizes_shelf) %>%mutate(community ="Wigley Subset",survey_area =factor(survey_area, levels = area_levels_all)) %>%left_join(area_df)# Get trends:wt_mod_wig <-gls( med_wt_kg ~scale(est_year) * area * season,correlation =corAR1(form =~ est_year | area/season),data =drop_na(size_results, med_wt_kg))# Get the predictions and flag whether they are significantsize_trend_predictions <-get_preds_and_trendsignif(wt_mod_wig) %>%mutate(model_id ="med_wt_kg-Wigley Subset") %>%separate(model_id, into =c("metric", "community"), sep ="-") %>%mutate(metric =fct_rev(metric))# Median length - all speciessize_long <- size_results %>%pivot_longer(cols =c(med_len_cm, med_wt_kg), names_to ="metric", values_to ="val") %>%mutate(survey_area =fct_relevel(survey_area, area_levels_all),area =fct_relevel(area, area_levels_long_all),season =fct_rev(season))# Median weight - wigley species# weight plotswt_plot <- size_long %>%filter(metric =="med_wt_kg") %>%drop_na(val) %>%ggplot() +geom_ribbon(data =filter(size_trend_predictions, metric =="med_wt_kg", non_zero ==TRUE),aes(x, ymin = conf.low, ymax = conf.high, fill = season),linewidth =1, alpha =0.35) +geom_point(aes(est_year, val, color = season),size =0.8, alpha =0.5) +geom_ma(aes(est_year, val, color = season, linetype ="5-Year Moving Average"), n =5, ma_fun = SMA, key_glyph ="timeseries") +geom_line(data =filter(size_trend_predictions, metric =="med_wt_kg", non_zero ==TRUE),aes(x, predicted, color = season, linetype ="Significant Trend"),linewidth =1, alpha =1) +scale_fill_gmri() +scale_color_gmri() +facet_grid(area ~ .,scales ="free") +labs(x ="Year",y ="Median Weight (kg)",color ="Season",linetype ="Trend",fill ="Season")wt_plot
Notes:
There are a number of years when median weight surges nearly 0.5kg in the Mid-Atlantic Bight. For these situations I think we likely could trace these down to exceptionally high catches of larger species like elasmobranchs.
We are not seeing a see a shelf-wide decline in size so I might need to check what Friedland did or prepare for a confrontation.
Bottom temperature data is from the DuPontavice and Saba combined bottom temperature product. This data product has been clipped to the study areas and seasonally averaged to match the survey sampling seasons.
Plotted below are the long-term trends in each area for each survey season. Significant seasonal trends are shown with linear trends and error displayed.
Seasonal bottom temperatures have been increasing in all 8/10 seasonal communities.
Spring temperature in SNE + MAB have not increased over this period.
Seasonal bottom temperatures have increased in both spring and fall for the Northeast Shelf.
Looking at different areas within the region, Fall temperatures are increasing across the board, but seasonally spring temperatures in SNE and MAB are stable.
NOTE: Percentile location of the size distribution are determined using DescTools::Quantile() which uses numlen_adj to weight everything by abundance, properly accounting for abundance.
The following figure shows changes in seasonal large and small fish indices within each region. The large and small size thresholds are set for each species independently using their length distributions as in the table above.
The axes on these are: \[LFI = \frac{\sum({N_{tjsk}
> size~threshold_k})}
{\sum{N_{tjs}}}\]
for each: \(year_t\), \(region_j\), \(season_s\) & \(species_k\).
Code
# Get percentiles based on length for each species# # Get 90th or 95th percentile as size threshold# Do them separatelyspecies_size_quants <- wigley_in %>%mutate(area ="Northeast Shelf") %>%unite("species_area", c(area, comname), sep ="X", ) %>%split(.$species_area) %>%map_dfr(function(x,y){ DescTools::Quantile(# x$length_cm, x$ind_weight_kg,weights = x$numlen_adj, probs =c(0.05, 0.1, 0.5, 0.90, .95)) %>%t() %>%as_tibble() }, .id ="species_area") %>%separate(col = species_area, into =c("area", "comname"), sep ="X") %>%mutate(area =factor(area, levels = area_levels_long_all)) %>%arrange(area)# # Show the table# species_size_quants %>% # mutate_if(is.numeric, round, 2) %>% # head(n = 10) %>% # gt() %>% # tab_header(title = "Wigley Community - Body Mass Quantiles (kg)") %>% tab_spanner(label = "Percentile", columns = -area) %>% # cols_label("area" = "")
# Plot them as loing-term average# Plot the SFIsfi_p <- lfi %>%ggplot() +geom_point(aes(est_year, SFI, color = season),size =0.8, alpha =0.5) +geom_ma(aes(est_year, SFI, color = season, linetype ="5-Year Moving Average"), n =5, ma_fun = SMA, key_glyph ="timeseries") +scale_fill_gmri() +scale_color_gmri() +scale_y_continuous(labels =label_percent()) +theme(strip.text.y =element_blank()) +facet_grid(area~., scales ="fixed", labeller =label_wrap_gen(width =8)) +labs(x ="Year",y ="Percent of Total Abundance",title ="Wigley Species",subtitle ="Small Fish Index",fill ="Season",color ="Season",linetype ="")# Plot the LFIlfi_p <-ggplot(lfi) +geom_point(aes(est_year, LFI, color = season),size =0.8, alpha =0.5) +geom_ma(aes(est_year, LFI, color = season, linetype ="5-Year Moving Average"), n =5, ma_fun = SMA, key_glyph ="timeseries") +scale_fill_gmri() +scale_color_gmri() +scale_y_continuous(labels =label_percent(), expand =expansion(add =c(0,0.05))) +facet_grid(area~., scales ="fixed", labeller =label_wrap_gen(width =8)) +labs(x ="Year",y ="",subtitle ="Large Fish Index",fill ="Season",color ="Season",linetype ="")# Plot together(sfi_p | lfi_p) +plot_layout(guides ="collect")
From these figures we can see that there is a general pattern where larger individuals (relative to their conspecific size-distribution) make up a smaller fraction of seasonal abundances.
There is something strange going on with GB in the 2000’s where the LFI increases for a time, but this is no longer the case.
Upper Percentile Drift of the Size distribution
These metrics are an attempt to capture where movements at the edges of the size distribution: Is the size distribution becoming wider or more narrow?
There is a sentiment that the 90th percentile fish is smaller now than it has been in the past. This is true in some circumstances but not all. Indeed, the 90th percentile mass for fishes in the GOM was nearly 2x what it is today. A similar change in the upper percentile location is visible in Georges Bank as well. At the shelf-wide scale and in the southern regions this is less true.
10th percentile weights seem stable and are likely less informative due to gear selectivity preventing consistent sampling at the smallest sizes. I find it more probable that the size distribution changes are related to fluctuations in larger individuals than those at the lower end of the distribution.
Code
# Do them separatelyyearly_size_quants <-bind_rows(mutate(wigley_in, area ="Northeast Shelf"), wigley_in) %>%unite(col ="group_var", area, est_year, sep ="X", remove = T) %>%split(.$group_var) %>%map_dfr(function(x,y){ DescTools::Quantile(# x$length_cm, x$ind_weight_kg,weights = x$numlen_adj, probs =c(0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.90, .95)) %>%t() %>%as_tibble() }, .id ="group_var") %>%separate(group_var, into =c("area", "est_year"), sep ="X") %>%mutate(est_year =as.numeric(est_year),area =factor(area, levels = area_levels_long_all)) %>%arrange(area)# Reshape to make filtering easierquants_long <- yearly_size_quants %>%pivot_longer(cols =-c(area, est_year), names_to ="Percentile", values_to ="biomass_kg") #-------------# Take the data excluding dogfishno_spiny <-filter(wigley_in, comname !="spiny dogfish")# Do them separatelyspineless_size_quants <-bind_rows(mutate(no_spiny, area ="Northeast Shelf"), no_spiny) %>%unite(col ="group_var", area, est_year, sep ="X", remove = T) %>%split(.$group_var) %>%map_dfr(function(x,y){ DescTools::Quantile( x$ind_weight_kg,weights = x$numlen_adj, probs =c(0.1, 0.5, 0.90)) %>%t() %>%as_tibble() }, .id ="group_var") %>%separate(group_var, into =c("area", "est_year"), sep ="X") %>%mutate(est_year =as.numeric(est_year),area =factor(area, levels = area_levels_long_all)) %>%arrange(area)# Reshape to make filtering easierspineless_quants_long <- spineless_size_quants %>%pivot_longer(cols =-c(area, est_year), names_to ="Percentile", values_to ="biomass_kg")
From exploratory analyses we know that spiny dogfish are a large fraction of the biomass in some seasonal communities. Because of this, their behavior has a large impact on the community size structure. However despite their outsized impact they are not the whole story and there is a danger in over-emphasizing their role as they are not strictly benthic and are known to migrate within and beyond the bounds of the survey sampling area.
By excluding spiny dogfish from the community the 2000’s era increase in 90th percentile size for the Gulf of Maine goes away. Possibly giving support to sentiments that spiny dogfish are occupying the body-size-niche perviously held by larger gadids.
The seasonality in the 90th percentile shift with and w/o spiny dogfish does not speak to migration as much as I thought it might originally.
It is capturing a more general absence of dogfish (and other larger individuals) in MAB during the fall. The departure of larger individuals from the area seasonally is possibly related to the breakdown of the Mid-Atlantic Cold pool, which breaks down over the summer.
Code
# Separate axes for 90quants_90 <- wigley_season_qs %>%filter(Percentile =="90%") %>%mutate(Percentile ="90th Percentile")spiny_90 <- spiny_season_qs %>%filter(Percentile =="90%") %>%mutate(Percentile ="90th Percentile")# 90th percentile dogfish changesq90_df_p <-ggplot() +geom_area(data = quants_90,aes(x = est_year, y = kg_roll, fill ="Complete Community"),alpha =0.6, color ="gray70") +geom_area(data = spiny_90,aes(x = est_year, y = kg_roll, fill ="Spiny Dogfish Removed"),alpha =0.6, color ="gray20") +geom_image(data =filter( quants_90, est_year ==2008, !(area =="Mid-Atlantic Bight"& season =="Fall")), aes(est_year, kg_roll*1.65, image = here::here("Data/phylopics/squalus_phylopic.png")),size =0.35) +scale_fill_manual(values =c("gray70", "gray20")) +scale_x_continuous(expand =expansion(add =c(0,0))) +scale_y_continuous(labels =label_number(suffix ="kg")) +#coord_flip() +guides(linetype =guide_legend(title.position ="top", title.hjust =0.5),color =guide_legend(title.position ="top",title.hjust =0.5,nrow =2)) +theme(plot.margin =margin(0,0,0,0),strip.background =element_rect(color ="white", linewidth =0.2)) +facet_nested( area~Percentile*season,labeller =labeller(area =label_wrap_gen(width =8),Percentile =~str_c(.x, " of Body Mass Distribution"))) +labs(y ="Individual Bodymass", x ="Year", fill ="Community")q90_df_p
The above figure shows that the community distribution without dogfish widens in the spring and narrows in the fall.
It is maybe worth noting just how much of the community biomass is held in spiny dogfish, which is frequently above 50%.
Spring bottom temperatures lead size distributions by 4 or 3 years. The spring temperature part may be important because that is the annual low i.e. the baseline is shifting.
Regionally:
Bottom temperature appears fairly centered (sequential significant correlations at or around 0 lag/lead) in GOM & GB where it is more highly correlated and there is more consistent correlation direction and magnitude at sequential lags.
For MAB-Spring a 1 year of bottom temperature may be important.
In the Fall the GOM has its largest negative correlation centered at a 2-year lead.
Code
# Function to grab the correlation data and lag dataget_ccf_vector <-function(x,y, lagmax =10){# Run the ccf ccf_data <-ccf(x, y, plot= F , lag.max = lagmax)# Get the signif:# https://www.squaregoldfish.co.uk/programming/r_acf_significance.md/# Not 100% sure n is the same for ccf as it is for acf, but... significance_level <-qnorm((1+0.95)/2)/sqrt(sum(!is.na(x)))data.frame("acf"= ccf_data$acf,"lag"= ccf_data$lag,"sigpos"= significance_level,"signeg"= significance_level*-1 )}# Run Local CCF on Bottom Temperaturebtemp_ccf <- wigley_bmspectra_model_df %>%unite(col ="regseas", area, season, sep ="XX") %>%split(.$regseas) %>%map_dfr(~get_ccf_vector(x = .x$bot_temp, y = .x$b), .id ="regseas") %>%separate(col ="regseas", into =c("area", "season"), sep ="XX") %>%mutate(# Flag when it crosses thresholdsig_flag =ifelse(acf < signeg | acf > sigpos, T, F),# Set Factor Levelsarea =factor(area, levels = area_levels_long_all),season =factor(season, levels =c("Spring", "Fall")))ggplot(btemp_ccf, aes(lag, acf)) +geom_hline(yintercept =0, color ="gray25", linewidth =1) +geom_vline(xintercept =0, color ="gray25", linewidth =1) +geom_ribbon(aes(ymin = signeg, ymax = sigpos), alpha =0.2, linetype =2, fill ="lightgray",color ="gray65") +geom_col(alpha =0.65) +geom_col(data =filter(btemp_ccf, sig_flag), color ="black") +scale_x_continuous(expand =expansion(add =c(0,0))) +facet_grid(area~season) +labs(title ="Bottom Temperature CCF",x ="Bottom Temperature Lag (years)", y ="ACF", )
Landings CCF
Its plausible that landings have a lagged effect so here are the CCFs for that.
These visually show a pattern of sustained autocorrelation, but in most instances the correlation to landings is such that the slope of the size distribution changes leads the landings AND that higher landings are positively correlated with shallow spectra.
Code
# Run Local CCF on Bottom Temperaturelandings_ccf <- wigley_bmspectra_model_df %>%drop_na(total_weight_lb) %>%unite(col ="regseas", area, season, sep ="XX") %>%split(.$regseas) %>%map_dfr(~get_ccf_vector(x = .x$total_weight_lb, y = .x$b), .id ="regseas") %>%separate(col ="regseas", into =c("area", "season"), sep ="XX") %>%mutate(# Flag when it crosses thresholdsig_flag =ifelse(acf < signeg | acf > sigpos, T, F),# Set Factor Levelsarea =factor(area, levels = area_levels_long_all),season =factor(season, levels =c("Spring", "Fall")))ggplot(landings_ccf, aes(lag, acf)) +geom_hline(yintercept =0, color ="gray25", linewidth =1) +geom_vline(xintercept =0, color ="gray25", linewidth =1) +geom_ribbon(aes(ymin = signeg, ymax = sigpos), alpha =0.2, linetype =2, fill ="lightgray",color ="gray65") +geom_col(alpha =0.65) +geom_col(data =filter(landings_ccf, sig_flag), color ="black") +scale_x_continuous(expand =expansion(add =c(0,0))) +facet_grid(area~season) +labs(title ="Total Landings CCF",x ="Total Landings Lag (years)", y ="ACF", )
Spectra Slope ACF
Just as an important check, what is the temporal autocorrelation structure in spectra slopes?
There is a significant positive correlation at 1 year lead in 4/10 regional timeseries. I think the nature of the data is enough to justify adding the AR1() error structure, and this is confirms its inclusion over a year random effect.
There is also six-year lead in 3 of them which is curious, and could be explored as to what is going on there with a 6-year periodicity.
Then lastly the Fall Gulf of Maine spectra has significant correlations up to three years out.
I am of two minds about this:
This could be a clue to some state change, as increased autocorrelation is sometimes used in the literature as a signal of hysteresis.
OR This could be evidence that Gulf of Maine’s fall community is the exception, and is the ONLY community sampled sufficiently/consistently.
We are in a sense (hopefully) conducting a repeated measure of some community in time, and so there should hopefully be some sort of autocorrelation. The absence of autocorrelation in these communities could
Code
# Pull out acf stuffspectra_acf <- wigley_bmspectra_model_df %>%unite("regseas", area, season, sep ="XX") %>%split(.$regseas) %>%map_dfr(function(x){# drop NA x <-arrange(x, est_year) %>%drop_na(b)# Do ACF x_acf <-acf(x$b, plot = F, lag.max =10)# Get the signif:# https://www.squaregoldfish.co.uk/programming/r_acf_significance.md/# Not 100% sure n is the same for ccf as it is for acf, but... significance_level <-qnorm((1+0.95)/2)/sqrt(sum(!is.na(x$b)))data.frame("acf"= x_acf$acf,"lag"= x_acf$lag,"sigpos"= significance_level,"signeg"= significance_level*-1 ) }, .id ="regseas") %>%separate("regseas", into =c("area", "season"), sep ="XX") %>%mutate(# Flag when it crosses thresholdsig_flag =ifelse(acf < signeg | acf > sigpos, T, F),# Set Factor Levelsarea =factor(area, levels = area_levels_long_all),season =factor(season, levels =c("Spring", "Fall"))) # Plot thingsggplot(spectra_acf, aes(lag, acf)) +geom_hline(yintercept =0, color ="gray25", linewidth =1) +geom_vline(xintercept =0, color ="gray25", linewidth =1) +geom_ribbon(aes(ymin = signeg, ymax = sigpos), alpha =0.2, linetype =2, fill ="lightgray",color ="gray65") +geom_col(alpha =0.65) +geom_col(data =filter(spectra_acf, sig_flag), color ="black") +scale_x_continuous(limits =c(-1,10),breaks =c(1:9),expand =expansion(add =c(0,0))) +facet_grid(area~season) +labs(title ="Spectra Slope ACF",x ="Lag (years)", y ="ACF", )
Spring and Fall Correlated ??
I can’t believe I never checked this, but are spring and fall communities correlated within an area? My sense is that they should be, and it might be meaningful to point out if they are not.
Only one region (GOM) has “memory” of the community six months before, and has a significant (but small) correlation between spring and fall communities.
Taken with the ACF info, there is consistency that the Gulf of Maine spectra exhibits more sign of temporal autocorrelation and the spectra between season are correlated.
This does however add doubt to how useful these other indices are, in the sense that they may not be capturing change of a community over time.
Size spectra steepening and median size declines are fairly consistent with one another. Long-term trends are only present at regional scales, and not ocurring shelf-wide for this community of 74 species.
There is something unique happening during the spring in the Mid-Atlantic Bight, where median size is growing and the size spectra is shallowing. This region exhibits large seasonal differences in community composition and large fluctuations in bottom temperature partly due to the mid-atlantic cold-pool.
From the large fish indices I see evidence that individuals within each species are smaller now than they had been, and that larger individuals are rarer in the overall community.
The shift in the 90th percentile size also supports this idea of smaller individuals being more common. The differences seen when excluding spiny dogfish gives some additional evidence that this is particularly true for non-dogfish species. This may also support the idea that dogfish species replacement has ocurred and spiny dogfish are now occupying that larger size class niche space.
From the temperature CCf plot my takeaways are that spring bottom temperature matters more for the shelf-wide size distribution and may operate with a lag of several years. Regionally, the lag seems less relevant and a negative correlation is more evident in GOM+GB.
From the landings CCF I get the sense that total landings aren’t a major driver of slope steepening in our area, or they operate at very long lead times 8-10 years. Many of the significant correlations are in the opposite direction of how we think landings would impact the size distribution (removal of larger individuals leads to steepening), AND they are such that they follow the size distribution slope rather than lead it.
Part 2: Relationships to External Factors
The aim of this section is to provide context and potentially attribution/correlations that may help explain trends/changes observed in the above section.
Temperature and Size Structure
There are known relationships between inter-specific and intra-specific body size across temperature ranges and across latitudinal ranges.
Bergmann’s Rule: Warmer climates are often inhabited by smaller-bodied species
James Rule: At an intraspecific level, populations living in warmer conditions often reach smaller maximum body sizes
The metabolic theory of ecology provides a basis for “why” the above rules might be observed.
Metabolic Theory of Ecology: Metabolic theory predicts how metabolic rate, by setting the rates of resource uptake from the environment and resource allocation to survival, growth, and reproduction, controls ecological processes at all levels of organization from individuals to the biosphere.
Seasonal Movements and Metabolic Efficiency: MTE details how metabolic rates, which are influenced by body size and temperature, are central to the ecological behavior of organisms.
This theory would suggest that larger individuals, which naturally have lower per-unit mass metabolic rates, may move to cooler areas seasonally to optimize their metabolic function and to reduce the energetic costs of thermoregulation in warmer periods.
From what I see at first-glance the size distribution along the shelf is responding consistently with temperature changes (Spring to Fall) on shorter time scales than hypotheses that predict temperature related changes in growth.
Code
wigley_bmspectra_model_df %>%filter(survey_area !="Northeast Shelf") %>%ggplot(aes(bot_temp, b, color = season)) +geom_point(aes(shape = survey_area)) +geom_smooth(method ="lm", color ="black") + ggforce::geom_mark_ellipse(aes(fill = season), alpha =0.1) +scale_color_gmri() +scale_fill_gmri() +scale_x_continuous(labels =label_number(suffix = deg_c)) +guides(shape =guide_legend(nrow =2),fill =guide_legend(nrow =2),color =guide_legend(nrow =2)) +theme(legend.box ="horizontal") +labs(title ="Bergmann's Rule + MTE",subtitle ="The broad relationship between temperature and size distribution is bolstered by\nseasonal community changes in alignment with seasonal temperature differences.",shape ="Survey Area",fill ="Season",color ="Season",x ="Bottom Temperature")
Importance of Seasonality
The extent that the regional spectra trends follow a negative linear relationship to increasing temperature is largely due to seasonal differences.
We know that individuals may migrate and/or occupy different habitats seasonally, leading to situations where the community compositions from spring to fall are an apples to oranges comparison. Seasonal movements may be be partially explained by expectations from the MTE, but they also could be separate ecological processes like seasonal foraging or spawning opportunities.
Either way, seasonal differences are large and they justify treating survey seasons as separate groups with independent trends.
Anomalies Model
I wanted to try and avoid conflating that a relationship between spectra and bottom temperature means that “warming” has altered the size spectra. The observation that body size distribution broadly follows some temperature relationship is a confirmation of these broadly observed rules, but I would argue it does not convey that warming has shifted a local climate in some direction.
This is the model I would suggest gets the most at “Warming” impacting regional size distributions:
Fixed effects:
Season
Survey area
Bottom temperature anomaly
Error structure: AR1 autocorrelative errors
By using the season * region interaction groups, we can characterize the typical bottom temperature for each combinations using the long-term seasonal average, and quantify warming using departures from those long-term averages.
The model looks like this as an equation (90% sure), where \(\alpha\) is the intercept, \(\beta x_{tjs}\) is our coefficient for a continuous independent variable (ex. bottom temperature), and \(e_{tjs\)` is our error:
\[y_{tjs} = \alpha + \beta x_{tjs} + e_{tjs}\] and our error depends on the previous timestep:
\[e_{tjs} = \phi e_{t-1,js} + w_{tjs}\] With \(w_{tjs}\) adding noise.
\[w_{tjs} \sim N(0, \sigma^2) \] All of this happening for each \(year_t\) within \(region_j\) & \(season_s\).
I need to consult the Zuur paper about equations in ecology again, but this seems close.
Code
# Just drop the random effect part for nowanom_model_ar <- nlme::gls( b ~ area * season *scale(bot_temp_anom), data = wigley_bmspectra_model_df, correlation =corAR1(form =~ est_year | area/season))# # confidence intervals on phi# intervals(anom_model_ar)$corStruct# # check the model# check_model(anom_model_ar)# plot(check_collinearity(anom_model_ar))# plot(# check_collinearity(# nlme::gls(# b ~ area + season + scale(bot_temp_anom), #+ scale(log(total_weight_lb)), # data = wigley_bmspectra_model_df, # correlation = corAR1(form = ~ est_year | area/season)# ))# )# Clean up the plot:# Plot the predictions over datamod2_preds <-as.data.frame(ggpredict( anom_model_ar, terms =~ bot_temp_anom * area * season)) %>%mutate(area =factor(group, levels = area_levels_long_all),season =factor(facet, levels =c("Spring", "Fall")))#### Trend Posthoc ##### Flag effect fits that are non-significant ###mod2_emtrend <-emtrends(object = anom_model_ar,specs =~ area * season,var ="bot_temp_anom") %>%as_tibble() %>%mutate(zero =0,non_zero =if_else(between(zero, lower.CL, upper.CL), F, T))# Limit temperature plotting rangeactual_range <- wigley_bmspectra_model_df %>%group_by(season, area) %>%summarise(min_temp =min(bot_temp_anom)-2,max_temp =max(bot_temp_anom)+2,.groups ="drop") %>%mutate(season =factor(season, levels =c("Spring", "Fall")))# Plot over observed datalocal_btemp_anoms <- mod2_preds %>%left_join(actual_range) %>%filter((x < min_temp) == F, (x > max_temp) == F) %>%left_join(select(mod2_emtrend, area, season, non_zero)) %>%filter(non_zero) %>%ggplot() +geom_ribbon(aes(x, ymin = conf.low, ymax = conf.high, group = season), alpha =0.1) +geom_line(aes(x, predicted, color = season), linewidth =1) +geom_point(data = wigley_bmspectra_model_df,aes(bot_temp_anom, b, color = season),alpha =0.4,size =1) +facet_grid(area~., scales ="free", labeller =label_wrap_gen(width =8)) +scale_color_gmri() +scale_x_continuous(labels =label_number(suffix = deg_c)) +labs(y ="Size Distribution Exponent",x ="Local Seasonal Bottom Temperature Anomaly")local_btemp_anoms
Based on this model I would say there is evidence of a warming effect in:
1. GOM-Spring
2. GB-Spring
3. GB-Fall
4. SNE-Fall
Bottom Temperature Model
We were originally using absolute temperatures, which are scaled, and I was curious with the AR error structure how different this would be than the anomaly model.
I just had to see if the same model is better performing or leads to different results when it knows absolute temperatures.
Code
# Just drop the random effect part for nowtemp_model_ar <- nlme::gls( b ~ area * season *scale(bot_temp), data = wigley_bmspectra_model_df, correlation =corAR1(form =~ est_year | area/season))# # confidence intervals on phi# intervals(temp_model_ar)$corStruct# # # check the model# check_model(temp_model_ar)# # Collinearity without interactions# plot(# check_collinearity(# nlme::gls(# b ~ area + season + scale(bot_temp), # data = wigley_bmspectra_model_df, # correlation = corAR1(form = ~ est_year | area/season))# ))# Clean up the plot:# Plot the predictions over datamod2_preds <-as.data.frame(ggpredict( temp_model_ar, terms =~ bot_temp * area * season)) %>%mutate(area =factor(group, levels = area_levels_long_all),season =factor(facet, levels =c("Spring", "Fall")))#### Trend Posthoc ##### Flag effect fits that are non-significant ###mod2_emtrend <-emtrends(object = temp_model_ar,specs =~ area * season,var ="bot_temp") %>%as_tibble() %>%mutate(zero =0,non_zero =if_else(between(zero, lower.CL, upper.CL), F, T))# Limit temperature plotting rangeactual_range <- wigley_bmspectra_model_df %>%group_by(season, area) %>%summarise(min_temp =min(bot_temp)-2,max_temp =max(bot_temp)+2,.groups ="drop") %>%mutate(season =factor(season, levels =c("Spring", "Fall")))# Plot over observed datalocal_btemp_results <- mod2_preds %>%left_join(actual_range) %>%filter((x < min_temp) == F, (x > max_temp) == F) %>%left_join(select(mod2_emtrend, area, season, non_zero)) %>%filter(non_zero) %>%ggplot() +geom_ribbon(aes(x, ymin = conf.low, ymax = conf.high, group = season), alpha =0.1) +geom_line(aes(x, predicted, color = season), linewidth =1) +geom_point(data = wigley_bmspectra_model_df,aes(bot_temp, b, color = season),alpha =0.4,size =1) +facet_grid(area~., scales ="free", labeller =label_wrap_gen(width =8)) +scale_color_gmri() +scale_x_continuous(labels =label_number(suffix = deg_c)) +labs(y ="Size Distribution Exponent",x ="Seasonal Bottom Temperature")local_btemp_results
Verdict: Temperature itself performs better, and is simpler to produce/explain.
Consistent seasonal shifts in the Shelf-wide spectra suggests to me that the community is more responsive to within-year time scale factors and maybe less sensistive to long-term trends in temperature than I originally expected. This also raises questions about how the spring/fall size distribution differences are being achieved and whether to view them as expressions of the same “community” or separately.
From the decadal variability work we have reason to believe that individual species are both not (perfectly) tracking temperature changes in space at the rate they are ocurring AND that changes in size at age are happening. These are signals at the species level indicating that temperature effects on growth are ocurring.
In light of that knowledge, there must be some compensatory mechanisms happening across this broader community for the temperature+spectra relationship to be so limited/situational.
Region * Season Pairwise Posthoc
The following figures and tables compare marginal means and post-hoc comparisons for the model above using bottom temperature.
What are the differences across spatial scales and seasonally?
Code
# Region and seasonal posthocsregseas_phoc <-emmeans( temp_model_ar, list(pairwise ~ area + season), adjust ="tukey",type ="response")# Custom Plot(intercepts_emmeans <- regseas_phoc$`emmeans of area, season`%>%as_tibble() %>%ggplot() +geom_pointrange(aes(area, emmean, ymin = lower.CL, ymax = upper.CL, color = season),position =position_dodge(width =0.25), alpha =1) +#coord_flip() +scale_color_gmri() +labs(y ="Marginal Mean: Mass Spectra Slope",x ="Area",title ="Wigley Species",color ="Season",subtitle ="Group marginal means plot"))
Code
# Table for pairwise test resultsregseas_phoc %>%pluck(2) %>%as.data.frame() %>%rename(Contrast =`1`) %>%mutate_if(is.numeric, ~round(.x, 3)) %>%arrange(p.value) %>%gt() %>%tab_header(title ="Factor Groups Pairwise Post-Hoc") %>%tab_style(style =list(cell_fill(color ="#9AFF9A"), cell_text(color ="black")),locations =cells_body(rows = p.value <=0.05) )
Factor Groups Pairwise Post-Hoc
Contrast
estimate
SE
df
t.ratio
p.value
Northeast Shelf Spring - (Mid-Atlantic Bight Fall)
0.943
0.130
471.630
7.252
0.000
Gulf of Maine Spring - (Mid-Atlantic Bight Fall)
0.878
0.134
475.331
6.549
0.000
Georges Bank Spring - (Mid-Atlantic Bight Fall)
0.862
0.150
474.065
5.731
0.000
Southern New England Spring - (Mid-Atlantic Bight Fall)
Southern New England Fall - (Mid-Atlantic Bight Fall)
0.742
0.125
476.724
5.913
0.000
(Mid-Atlantic Bight Spring) - Northeast Shelf Fall
0.245
0.061
458.918
4.011
0.003
Georges Bank Fall - Southern New England Fall
0.375
0.094
475.495
3.976
0.003
Northeast Shelf Fall - Gulf of Maine Fall
-0.292
0.077
460.546
-3.815
0.006
Gulf of Maine Fall - Southern New England Fall
0.283
0.089
471.001
3.173
0.051
(Mid-Atlantic Bight Spring) - Southern New England Fall
0.237
0.077
471.559
3.091
0.065
Southern New England Spring - Northeast Shelf Fall
0.198
0.079
470.459
2.518
0.262
Gulf of Maine Spring - Georges Bank Fall
-0.239
0.106
477.157
-2.261
0.417
Northeast Shelf Spring - Northeast Shelf Fall
0.209
0.094
463.762
2.226
0.441
Southern New England Spring - Georges Bank Fall
-0.186
0.086
477.144
-2.156
0.489
Southern New England Spring - Southern New England Fall
0.190
0.091
474.120
2.079
0.543
Georges Bank Spring - Georges Bank Fall
-0.255
0.126
476.856
-2.031
0.578
(Mid-Atlantic Bight Spring) - Georges Bank Fall
-0.139
0.070
476.051
-1.970
0.621
Northeast Shelf Spring - Southern New England Fall
0.201
0.105
469.083
1.917
0.657
Northeast Shelf Spring - Georges Bank Fall
-0.174
0.100
471.737
-1.739
0.773
Gulf of Maine Spring - Northeast Shelf Fall
0.145
0.100
473.576
1.454
0.909
Gulf of Maine Spring - Gulf of Maine Fall
-0.147
0.101
476.177
-1.453
0.909
Georges Bank Spring - Gulf of Maine Fall
-0.163
0.122
476.375
-1.340
0.944
Gulf of Maine Spring - Southern New England Fall
0.136
0.110
471.249
1.243
0.965
Southern New England Spring - Gulf of Maine Fall
-0.094
0.081
477.748
-1.165
0.977
Gulf of Maine Spring - (Mid-Atlantic Bight Spring)
-0.100
0.090
473.850
-1.112
0.983
Gulf of Maine Fall - Georges Bank Fall
-0.092
0.084
444.940
-1.092
0.985
Georges Bank Spring - Northeast Shelf Fall
0.129
0.121
469.365
1.066
0.988
Georges Bank Spring - (Mid-Atlantic Bight Spring)
-0.116
0.113
477.514
-1.032
0.990
Georges Bank Spring - Southern New England Fall
0.120
0.129
431.576
0.930
0.995
Northeast Shelf Spring - Gulf of Maine Fall
-0.083
0.096
469.882
-0.864
0.997
Southern New England Spring - (Mid-Atlantic Bight Spring)
-0.047
0.066
477.952
-0.711
0.999
(Mid-Atlantic Bight Spring) - Gulf of Maine Fall
-0.047
0.063
470.330
-0.738
0.999
Northeast Shelf Spring - Gulf of Maine Spring
0.064
0.115
477.403
0.560
1.000
Northeast Shelf Spring - Georges Bank Spring
0.081
0.134
475.084
0.604
1.000
Northeast Shelf Spring - Southern New England Spring
0.011
0.097
408.990
0.115
1.000
Northeast Shelf Spring - (Mid-Atlantic Bight Spring)
-0.036
0.084
474.661
-0.427
1.000
Gulf of Maine Spring - Georges Bank Spring
0.016
0.138
433.350
0.119
1.000
Gulf of Maine Spring - Southern New England Spring
-0.053
0.103
477.368
-0.517
1.000
Georges Bank Spring - Southern New England Spring
-0.070
0.123
476.681
-0.564
1.000
Northeast Shelf Fall - Southern New England Fall
-0.008
0.088
417.718
-0.096
1.000
Code
regseas_phoc %>%pluck(2) %>%as.data.frame() %>%rename(Contrast =`1`) %>%mutate_if(is.numeric, ~round(.x, 3)) %>%filter(p.value <0.05) %>%separate(Contrast, into =c("group1", "group2"), sep =" - ") %>%mutate(group1 =str_remove_all(group1, "[(]|[)]"),group2 =str_remove_all(group2, "[(]|[)]"),statement =if_else(estimate >0, " is shallower than ", " is steeper than ")) %>%select(group1, statement, group2) %>%arrange(statement, group1) %>%gt()
group1
statement
group2
Georges Bank Fall
is shallower than
Southern New England Fall
Georges Bank Fall
is shallower than
Mid-Atlantic Bight Fall
Georges Bank Spring
is shallower than
Mid-Atlantic Bight Fall
Gulf of Maine Fall
is shallower than
Mid-Atlantic Bight Fall
Gulf of Maine Spring
is shallower than
Mid-Atlantic Bight Fall
Mid-Atlantic Bight Spring
is shallower than
Northeast Shelf Fall
Mid-Atlantic Bight Spring
is shallower than
Mid-Atlantic Bight Fall
Northeast Shelf Fall
is shallower than
Mid-Atlantic Bight Fall
Northeast Shelf Spring
is shallower than
Mid-Atlantic Bight Fall
Southern New England Fall
is shallower than
Mid-Atlantic Bight Fall
Southern New England Spring
is shallower than
Mid-Atlantic Bight Fall
Northeast Shelf Fall
is steeper than
Gulf of Maine Fall
Northeast Shelf Fall
is steeper than
Georges Bank Fall
Total Landings
The landings information is the newer GARFO data obtained from Andy Beet. It was delivered as containing only finfish, I have removed freshwater species and blue crabs. It should largely reflect marine finfish species at this point.
I wanted to follow on the bottom temperature anomalies model here for landings:
Code
# Just drop the random effect part for nowf_model_ar <- nlme::gls( b ~ area * season + area *scale(log(total_weight_lb)), data =drop_na(wigley_bmspectra_model_df, total_weight_lb) %>% dplyr::filter(est_year >1978),correlation =corAR1(form =~ est_year | area/season))# Clean up the plot:# Plot the predictions over datamod2_preds <-as.data.frame(ggpredict( f_model_ar, terms =~ total_weight_lb * area * season)) %>%mutate(area =factor(group, levels = area_levels_long_all),season =factor(facet, levels =c("Spring", "Fall")))#### Trend Posthoc ####trend_df <-emtrends(object = f_model_ar,~area,var ="total_weight_lb",mode ="appx-satterthwaite")# Flag effect fits that are non-significant ###mod2_emtrend <-emtrends(object = f_model_ar,specs =~ area,var ="total_weight_lb",mode ="appx-satterthwaite") %>%as_tibble() %>%mutate(zero =0,non_zero =if_else(between(zero, lower.CL, upper.CL), F, T))# Limit temperature plotting rangeactual_range <- wigley_bmspectra_model_df %>%group_by(season, area) %>%summarise(min_f =min(total_weight_lb)-2,max_f =max(total_weight_lb)+2,.groups ="drop") %>%mutate(season =factor(season, levels =c("Spring", "Fall")))# Plot over observed datalandings_ar_plot <- mod2_preds %>%left_join(actual_range) %>%filter((x < min_f) == F, (x > max_f) == F) %>%left_join(select(mod2_emtrend, area,non_zero)) %>%filter(non_zero) %>%ggplot() +geom_ribbon(aes(x, ymin = conf.low, ymax = conf.high, group = season), alpha =0.1) +geom_line(aes(x, predicted, color = season), linewidth =1) +geom_point(data = wigley_bmspectra_model_df,aes(total_weight_lb, b, color = season),alpha =0.4,size =1) +facet_grid(area~., scales ="free") +scale_color_gmri() +scale_x_log10(labels =label_log(base =10), limits =10^c(0,10)) +labs(y ="Exponent of ISD",title ="Total Landings (lb.)",x ="Local Seasonal Bottom Temperature Anomaly")landings_ar_plot
Code
# temp model is better# AIC(temp_model_ar)# AIC(f_model_ar)
This is our old friend “higher landings when spectra was shallower” that we’ve seen before.
Temperature and Landings Together
If we want to just throw them both in because we’ve done what we can and want to hold them both in the frame at once we get this:
Code
# Just drop the random effect part for nowfull_model_ar <- nlme::gls(# Modelmodel = b ~ area * season *scale(bot_temp) + area *scale(log(total_weight_lb)), # Datadata = wigley_bmspectra_model_df %>%drop_na(total_weight_lb) %>% dplyr::filter(area !="Northeast Shelf"),# Autocorrelationcorrelation =corAR1(form =~ est_year | area/season))# # Collinearity without interactions# plot(# check_collinearity(# nlme::gls(# b ~ area + season + scale(bot_temp) + scale(log_land),# data = drop_na(wigley_bmspectra_model_df, total_weight_lb) %>%# filter(area != "Northeast Shelf") %>% # mutate(log_land = log(total_weight_lb)),# correlation = corAR1(form = ~ est_year | survey_area/season)# )# )# )# # confidence intervals on phi# intervals(full_model_ar)$corStruct# temp model is much better# AIC(temp_model_ar)# AIC(full_model_ar)
Before proceeding it is worth noting that this model is less parsimonious than a temperature or anomaly only model:
There are several cold pool indices available that might be helpful in adding context to changes in MAB. Below are the Dupontavice et al. 2022 methods that are used in the SOE report.
The Cold Pool Index (Model_CPI) was adapted from Miller et al. (2016). Residual temperature was calculated in each grid cell, \(i\), in the Cold Pool domain as the difference between the average bottom temperature at the year \(y\) (\(T_{i,y}\)) and the average bottom temperature over the period 1972–2019 (\(\bar{T}_{i, 1972-2019}\)) between June and September. \(CPI\) was calculated as the mean residual temperature over the Cold Pool domain such that